pica1_after_pica2 <- function(Y, D, A, C, S, 
                              placebo, placebo_level, data, 
                              # effects.focus,
                              posterior = TRUE,
                              boot = FALSE,
                              num_boot = 100,
                              n_dir_MC = 100,
                              n_norm_MC = 10,
                              quantiles = c(.025,.975),
                              sensitivity = FALSE,
                              rhos = Inf,
                              hpd = TRUE,
                              alpha = .95,
                              seed = 1234){
  
  # Cleaning data
  # remove   <- apply(data[, c(Y, D, A, S, placebo)], 1, function(x) any(is.na(x)))
  # data_use <- data[remove == FALSE, ]  
  
  set.seed(seed)
  data_use <- data
  
  data_use$Y <- data_use[, Y]
  data_use$D <- data_use[, D]
  data_use$A <- data_use[, A]
  data_use$C <- data_use[, C]
  data_use$S <- data_use[, S]
  data_use$placebo <- data_use[, placebo]
  
  # Analyze as the usual PICA (One key change)
  data_use$Y[data_use[, placebo] == 1] <- NA   
  J <- length(unique(na.omit(data_use$C)))
  
  # ATE 
  effects.mat_use <- matrix(FALSE, J, J)
  dimnames(effects.mat_use) <- list(ctrl = NULL, treat = NULL)
  effects.mat_use[placebo_level, ] <- TRUE
  effects.mat_use[placebo_level, placebo_level] <- FALSE
  
  # ACTE (Additional ACTE we have to merge to estimate)
  effects.focus <- cbind("treat" = seq(1:J)[-placebo_level], "ctrl" = rep(placebo_level, J -1), "C" = - seq(1:J)[-placebo_level])
  
  out_pica1 <- bounds.frechet_pica1_after_pica2(data = data_use, 
                                                J = J, 
                                                effects.mat = effects.mat_use,
                                                effects.focus = effects.focus, 
                                                posterior = posterior,
                                                n_dir_MC = n_dir_MC,
                                                n_norm_MC = n_norm_MC,
                                                quantiles = quantiles,
                                                sensitivity = sensitivity,
                                                rhos = rhos,
                                                hpd = hpd,
                                                alpha = alpha,
                                                verbose = 0  ## currently has no effect
  )
  
  if(posterior == TRUE){
    effects.ci_use <- out_pica1$posterior$effects.ci
    effects_out <- cbind(out_pica1$effects, effects.ci_use[, -c(1,2,3)])
    effects_out <- rbind(effects_out, out_pica1$effects.new_tab)
    
  }else if(boot == TRUE){
    out_pica1_boot <- bounds.frechet.boot_pica12(data_use_boot = data_use, num_boot = num_boot)
    effects.ci_use <- out_pica1_boot$effects.ci
    effect_new <- cbind(cbind("C" = c(-1, -2), "treat" = c(1, 2), "ctrl" = c(3,3)), out_pica1_boot$effect_new)
    effects_out <- rbind(cbind(out_pica1$effects, effects.ci_use[, -c(1,2,3)]), effect_new)
  }
  colnames(effects_out)[4:7] <- c("Lower Bound", "Upper Bound", "Lower CI", "Upper CI")
  rownames(effects_out) <- NULL
  return(effects_out)
}




bounds.frechet.boot_pica12 <- function(data_use_boot, num_boot = 100){
  
  ## bootstrapping
  bounds.boot.list = lapply(1:num_boot, function(j){
    # cat('bounds boot', j, '\n')
    tryCatch({
      bounds.frechet_pica1_after_pica2(data_use_boot[sample.int(nrow(data_use_boot), replace = TRUE),],
                                       J = J, 
                                       effects.mat = effects.mat,
                                       posterior=FALSE,
                                       sensitivity = TRUE,
                                       rhos = seq(0, 12, 1)
      )
    },
    error = function(e){
      NA
    })
  })
  
  ## process bootstrap results and convert to ci
  
  failed_bootstrap_draw = sapply(bounds.boot.list, length) == 1
  bounds.boot.list = bounds.boot.list[!failed_bootstrap_draw]
  
  ## naive strata
  naive_strata.boot.arr = sapply(
    bounds.boot.list,
    function(x){
      as.matrix(x$naive[['strata']])
    },
    simplify = 'array'
  )
  naive_strata.boot <- as.data.frame(apply(naive_strata.boot.arr, 1:2, quantile, .025))
  colnames(naive_strata.boot)[match('naive', colnames(naive_strata.boot))] <-
    'min_cilo'
  naive_strata.boot$max_cihi <- apply(naive_strata.boot.arr[,'naive',], 1, quantile, .975)
  
  ## naive effects
  naive_effects.boot.arr = sapply(
    bounds.boot.list,
    function(x){
      as.matrix(x$naive[['effects']])
    },
    simplify = 'array'
  )
  naive_effects.boot <- as.data.frame(apply(naive_effects.boot.arr, 1:2, quantile, .025))
  colnames(naive_effects.boot)[match('naive', colnames(naive_effects.boot))] <-
    'min_cilo'
  naive_effects.boot$max_cihi <- apply(naive_effects.boot.arr[,'naive',], 1, quantile, .975)
  
  ## bounds on strata
  strata.boot.arr = sapply(
    bounds.boot.list,
    function(x){
      as.matrix(x[['strata']])
    },
    simplify = 'array'
  )
  strata.boot <- as.data.frame(apply(strata.boot.arr, 1:2, quantile, .025))
  strata.boot$max <- apply(strata.boot.arr[,'max',], 1, quantile, .975)
  colnames(strata.boot)[match(c('min', 'max'), colnames(strata.boot))] <-
    c('min_cilo', 'max_cihi')
  
  ## bounds on effects
  effects.boot.arr = sapply(
    bounds.boot.list,
    function(x){
      as.matrix(x[['effects']])
    },
    simplify = 'array'
  )
  effects.boot <- as.data.frame(apply(effects.boot.arr, 1:2, quantile, .025))
  effects.boot$max <- apply(effects.boot.arr[,'max',], 1, quantile, .975)
  colnames(effects.boot)[match(c('min', 'max'), colnames(effects.boot))] <-
    c('min_cilo', 'max_cihi')
  
  ## NEW Estimate 
  bounds.boot.new_13_m1 <- bounds.boot.new_23_m2 <- matrix(NA, nrow = length(bounds.boot.list), ncol = 2)
  for(x in 1:length(bounds.boot.list)){
    # cat(x)
    if(any(is.na(bounds.boot.list[[x]])) == FALSE){
      
      prop_c1 <- bounds.boot.list[[x]]$prop_c[-1]/sum(bounds.boot.list[[x]]$prop_c[-1])
      prop_c2 <- bounds.boot.list[[x]]$prop_c[-2]/sum(bounds.boot.list[[x]]$prop_c[-2])
      
      bounds.boot.new_13_m1[x, 1:2] <- c(t(bounds.boot.list[[x]]$effects[c(5, 6), -c(1,2,3), x]) %*% 
                                           prop_c1)
      bounds.boot.new_23_m2[x, 1:2] <- c(t(bounds.boot.list[[x]]$effects[c(7, 9), -c(1,2,3), x]) %*% 
                                           prop_c2)
    }
  }
  effect_13_m1 <- c(apply(bounds.boot.new_13_m1[,c(1,2)], 2, mean, na.rm = TRUE),
                    quantile(bounds.boot.new_13_m1[,1], prob = 0.025, na.rm = TRUE),
                    quantile(bounds.boot.new_13_m1[,2], prob = 0.975, na.rm = TRUE))
  effect_23_m2 <- c(apply(bounds.boot.new_23_m2[,c(1,2)], 2, mean, na.rm = TRUE),
                    quantile(bounds.boot.new_23_m2[,1], prob = 0.025, na.rm = TRUE),
                    quantile(bounds.boot.new_23_m2[,2], prob = 0.975, na.rm = TRUE))
  effect_new <- rbind(effect_13_m1, effect_23_m2)
  rownames(effect_new) <- c("13_m1", "23_m2")
  colnames(effect_new) <- c("min", "max", "min_cilo", "max_cihi")
  
  ## sensitivity on strata
  sens_strata.boot.arr = sapply(
    bounds.boot.list,
    function(x){
      as.matrix(x$sensitivity[['strata']])
    },
    simplify = 'array'
  )
  sens_strata.boot <-
    as.data.frame(apply(sens_strata.boot.arr, 1:2, quantile, .025, na.rm = TRUE))
  sens_strata.boot$max <-
    apply(sens_strata.boot.arr[,'max',], 1, quantile, .975, na.rm = TRUE)
  colnames(sens_strata.boot)[match(c('min', 'max'), colnames(sens_strata.boot))] <-
    c('min_cilo', 'max_cihi')
  sens_strata.boot$failed_lp <-
    apply(sens_strata.boot.arr[,'min',], 1, function(x) mean(is.na(x)))
  
  ## sensitivity on effects
  sens_effects.boot.arr = sapply(
    bounds.boot.list,
    function(x){
      as.matrix(x$sensitivity[['effects']])
    },
    simplify = 'array'
  )
  sens_effects.boot <-
    as.data.frame(apply(sens_effects.boot.arr, 1:2, quantile, .025, na.rm = TRUE))
  sens_effects.boot$max <-
    apply(sens_effects.boot.arr[,'max',], 1, quantile, .975, na.rm = TRUE)
  colnames(sens_effects.boot)[match(c('min', 'max'), colnames(sens_effects.boot))] <-
    c('min_cilo', 'max_cihi')
  sens_effects.boot$failed_lp <-
    apply(sens_effects.boot.arr[,'min',], 1, function(x) mean(is.na(x)))
  
  return(list(
    naive_strata.ci = naive_strata.boot,
    naive_effects.ci = naive_effects.boot,
    strata.ci = strata.boot,
    effects.ci = effects.boot,
    sens_strata.ci = sens_strata.boot,
    sens_effects.ci = sens_effects.boot,
    failed_boot = mean(failed_bootstrap_draw),
    effect_new = effect_new
  ))
  
}
